% this is the script for calculating the fractal dimention (D-value) in 3D
% Version: 1.1
clc;
clear all;
close all;

%Read circle data file
f_id=fopen('2D-section.dat');
XYR=fscanf(f_id,'%f %f %f\n',[3 inf]);
f_return=fclose(f_id);

%X=XYZR(1,:);
%Y=XYZR(2,:);
R=XYR(3,:);
D=R*2;

maxD=max(R)*2;
minD=min(R)*2;
%NumSphere=size(R,2);

%Calculate D-value

binInterval=1;

tempBin=[-10:binInterval:10];
maxDLog2=log2(maxD);
maxDLog2=min(find(tempBin>maxDLog2));
minDLog2=log2(minD);
minDLog2=max(find(tempBin<minDLog2));

phiBin=seqreverse([tempBin(minDLog2):binInterval:tempBin(maxDLog2)]);
numPhiBin=size(phiBin,2)-1;
Dref=1;
DBin=Dref*power(2,phiBin);

cumulCount=0;

for i=1:numPhiBin
    iBin=find(D<=DBin(i) & D>DBin(i+1));
    count=size(iBin,2);
    cumulCount=cumulCount+count;
    cumulFreq(i)=cumulCount;
    if count==0
        diameter(i)=(DBin(i)+DBin(i+1))/2;
    else
        diameter(i)=median(D(iBin));
    end;
end;

rangeInit=1;
rangeFinal=numPhiBin;
rangeFinal=6;

logCFreq=log10(cumulFreq(rangeInit:rangeFinal));
logDiameter=log10(diameter(rangeInit:rangeFinal));
%plot(logDiameter,logCFreq,'+');
bestFit=polyfit(logDiameter,logCFreq,1);
D_Value=bestFit(1)*-1

%Plot the Log-Log diagram 
loglog(diameter(rangeInit:rangeFinal),cumulFreq(rangeInit:rangeFinal),'+');
axis_font_size=10;
xlabel('Diameter (log)','fontsize',axis_font_size);
ylabel('Cumulative Frequency (log)','fontsize',axis_font_size);


